clc
clear
currentDir = pwd; % 获取当前目录路径
parentDir = fileparts(currentDir); % 获取当前目录的上一级目录路径
% grandParentDir = fileparts(parentDir); % 获取当前目录的上两级目录路径
% greatGrandParentDir = fileparts(grandParentDir); % 获取当前目录的上三级目录路径

%指定需要查找的文件夹
subdir = 'DATA';
%指定需要查找的文件名
filename = '413168000.mat';
%构造需要查找的文件的完整路径
full_path = fullfile(parentDir, subdir, filename);

load( full_path)  %此处使用的数据是原始org数据经过后面的[]处理后得到的数据

%% 数据预处理
AIS(AIS.Speed<1,:)=[];%锚泊状态不处理

%% 重复值剔除
index_repeat=delete_repeat_value(AIS);
AIS(index_repeat,:)=[];


current_ais = AIS;
time_error_org={};%需要插值点数据time值
count_1=size(current_ais,1);
mean_pz=[];%保存前后窗口+异常数据集的合并数据集各阶段的平均偏转角，第一列为原始数据，第二列为剔除异常值后的数据，第三列为加入插值后的数据
data_interp=[];
while(1)  %如果连续多个点是异常数据，每次只处理一个，因此最好多次预处理
    [ais_deleted,time_error_org{end+1},current_pz,current_data_interp] = preprocess_AIS(current_ais);%剔除异常值后的数据
    mean_pz=[mean_pz;current_pz];
    count_2 = size(ais_deleted,1);
    data_interp=[data_interp;current_data_interp];
    if count_1-count_2 <3%当前一次处理与后一次处理的数据相差不大，则停止处理
        break;
    else
        current_ais=ais_deleted;
    end

    count_1 = count_2;
end

spline_time=vertcat(time_error_org{:});%合并数值，需要插值点数据time值
%需要将插值中的0值剔除
spline_time(spline_time==0,:)=[];%真正需要进行插值点数据对应的time值
spline_time = sort(spline_time);
time_diff = ais_deleted.time-ais_deleted.time(1);%与第一个数据的时间差值
time_interp = spline_time-ais_deleted.time(1);%待插值数据与第一个数据的时间差值

spline_data = AIS(ismember(AIS.time,spline_time),:);%需要插值点的原来数据
spline_data.Lon=spline(time_diff,ais_deleted.Lon,time_interp);%用插值替换掉原来数据
spline_data.Lat=spline(time_diff,ais_deleted.Lat,time_interp);

%将插值数据放回原来位置，即将剔除异常点的数据与插值数据进行结合
ais_splined=vertcat(spline_data,ais_deleted);%插值后的轨迹数据
%按照时间time值进行排序
ais_splined=sortrows(ais_splined,"time");

%绘制偏转角对比曲线并计算偏转角标准差,以及计算轨迹曲率
draw_process_result(AIS,ais_deleted,ais_splined,mean_pz);

figure
degrees_test = calculate_track_degree([ais_splined.Lon,ais_splined.Lat]);
plot(degrees_test);

% %绘制基于窗口的插值结果,从结果看，基于窗口过程的插值，收到的影响较大，产生的误差也较大
% figure
% window_spline=vertcat(ais_deleted,data_interp);%插值数据是基于窗口数据集的
% window_spline=sortrows(window_spline,"time");
% window_spline_degrees=calculate_track_degree([window_spline.Lon,window_spline.Lat]);
% plot(window_spline_degrees);

%% 数据集划分
%将原始数据直接划分为训练集和测试集，然后在分别对其进行重构
count_ais=size(AIS,1);
idx_train=1:floor(0.9*count_ais);
idx_test=floor(0.9*count_ais)+1:count_ais;
data_train=table2array( AIS(idx_train,:));
data_test=table2array(AIS(idx_test,:));
count_row_train=size(data_train,1);%训练集的行数
count_row_test=size(data_test,1);%测试集的行数
%分别重构训练集和测试集，时间步长为5
time_step=5;
train_x=cell(count_row_train-time_step,1);
train_y=zeros(size(train_x,1),count_col);

test_x=cell(count_row_test-time_step,1);
test_y=zeros(size(test_x,1),count_col);

%重构训练集
for i=1:count_row_train-time_step
    train_x{i}=data_train(i:i+time_step-1,:);
    train_y(i,:)=data_train(i+time_step,:);
end

%重构测试集
for i=1:count_row_test-time_step
    test_x{i}=data_test(i:i+time_step-1,:);
    test_y(i,:)=data_test(i+time_step,:);
end

%% 数据归一化,是每个维度数据的归一化
muX = mean(cell2mat(train_x));
sigmaX = std(cell2mat(train_x),0);

muY = mean(train_y);
sigmaY = std(train_y,0);
%训练数据归一化
for n = 1:numel(train_x)
    train_x{n} = (train_x{n} - muX) ./ sigmaX;
    train_y(n,:) = (train_y(n,:) - muY) ./ sigmaY;
end

%测试数据归一化
for n=1:numel(test_x)
    test_x{n} = (test_x{n} - muX) ./ sigmaX;
    %test_y(n,:) = (test_y(n,:) - muY) ./ sigmaY;
end

%% 定义LSTM网络结构
num_Features =size(train_x{1},2);   %输入特征维度
num_Responses = 1;%size(train_y,2);  %输出特征维度   可以对比单维度输出和5个维度一起输出时的结果
layers = [
    sequenceInputLayer(num_Features)
    lstmLayer(128,OutputMode="last")
    fullyConnectedLayer(num_Responses)];

%% 指定训练选项
options = trainingOptions("adam", ...
    MaxEpochs=100, ...
    SequencePaddingDirection="left", ...
    Shuffle="every-epoch", ...
    Plots="training-progress", ...
    Metrics="rmse",...
    Verbose=false);
%TargetDataFormats="CB", ... % 根据目标数据的实际格式设置
%% 训练网络
% y1=cellfun(@(x) x(:, 1:5), train_y, 'UniformOutput', false);
% y1=cell2mat(y1);
% y2=cellfun(@(x) x(:, 2), train_y, 'UniformOutput', false);
% y2=cell2mat(y2);
% y=[y1 y2];

for i=1:2  %分别预测经纬度数值
    net = trainnet(train_x,train_y(:,i),layers,"mse",options);
    %% 测试网络
    YTest = minibatchpredict(net,test_x, ...
        SequencePaddingDirection="left", ...
        UniformOutput=false);
    data_predict(:,i)= cell2mat(YTest);
    net=resetState(net);
end

%% 将获得的预测经纬度进行还原
data_predict=data_predict.*sigmaX(1:2)+muX(1:2);

plot(data_predict(:,1),data_predict(:,2),'k-o');
plot(test_y(:,1),test_y(:,2),'r*');

%
% % 在直方图中可视化均方误差。
% figure
% TTest = cellfun(@(x) x(:, 1), test_y, 'UniformOutput', false);
% TTest=cell2mat(TTest);
% YTest=cell2mat(YTest);
% histogram(mean((TTest -YTest ).^2,2))
% xlabel("Error")
% ylabel("Frequency")
%
% % 计算总体均方根误差。
% rmse = rmse(YTest,TTest);
%
% % 绘制预测频率对实际频率的图。
% figure
% scatter(YTest,TTest, "b+");
% xlabel("Predicted Frequency")
% ylabel("Actual Frequency")
% hold on
%
% m = min(y1);
% M=max(y1);
% xlim([m M])
% ylim([m M])
% plot([m M], [m M], "r--")